Install necessary packages:
install.packages("tidyverse", dependencies = TRUE)
install.packages("raster", dependencies = TRUE)
install.packages("sf", dependencies = TRUE)
install.packages("landscapetools", dependencies = TRUE)
install.packages("landscapemetrics", dependencies = TRUE)
Import only the RGB color bands as individual
RasterLayer objects:
library(raster)
#blue
b2 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B2.tif')
#green
b3 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B3.tif')
#red
b4 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B4.tif')
Combine the RasterLayer objects and visualise the
satellite image:
landsatRGB <- stack(b4, b3, b2) #order is impt
plotRGB(landsatRGB,
stretch = "lin") #scale the values (try using "hist" also)
Import all 5 bands from the satellite data as a
RasterStack object named landsat:
filenames <- paste0('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B', 1:5, ".tif")
landsat <- stack(filenames)
#rename bands
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR')
Check coordinate reference system of landsat:
crs(landsat)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 48N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",105,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16048]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]
Import polygon of city boundaries as sgshp and check if
their coordinate reference systems match:
library(sf)
sgshp <- st_read("data/master-plan-2014-region-boundary-web-shp/MP14_REGION_WEB_PL.shp")
Check coordinate reference system of sgshp:
crs(sgshp)
## [1] "PROJCRS[\"SVY21\",\n BASEGEOGCRS[\"SVY21[WGS84]\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",6326]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n CONVERSION[\"unnamed\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",1.36666666666667,\n ANGLEUNIT[\"Degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",103.833333333333,\n ANGLEUNIT[\"Degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",1,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",28001.642,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",38744.572,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"
Transform sgshp to the match the coordinate reference
system of the landsat:
sgshp <- st_transform(sgshp, crs = crs(landsat))
Crop landsat according to city boundaries:
landsat <- crop(landsat, sgshp) #crop to rectangle
landsat <- mask(landsat, sgshp) #mask values according to shape of sgshp
Plot the cropped image using only the RGB bands:
landsatRGB <- subset(landsat, c(4,3,2)) #Red, Green, Blue
plotRGB(landsatRGB,
stretch = "lin")
Create a function that calcuates the Normalized Difference Vegetation Index (NDVI) for each pixel:
ndvi <- function(x, y) {
(x - y) / (x + y)
}
Apply function to the NIR and Red bands of landsat
landsatNDVI <- overlay(landsat[[5]], landsat[[4]],
fun = ndvi)
Limit the range of values to be from -1 to 1:
landsatNDVI <- reclassify(landsatNDVI, c(-Inf, -1, -1)) # <-1 becomes -1
landsatNDVI <- reclassify(landsatNDVI, c(1, Inf, 1)) # >1 becomes 1
Map out the NDVI values:
plot(landsatNDVI,
col = rev(terrain.colors(10)),
main = "Landsat 8 NDVI",
axes = FALSE, box = FALSE)
Plot histogram of NDVI values:
hist(landsatNDVI,
main = "Distribution of NDVI values", xlab = "NDVI",
xlim = c(-1, 1), breaks = 100, yaxt = 'n')
abline(v=0.2, col="red", lwd=2)
abline(v=0, col="red", lwd=2)
Set 0.2 as the threshold; reclassify values below this threshold to
NA:
landsatGreen <- reclassify(landsatNDVI, c(-1, 0.2, NA)) #-1 to 0.2 becomes NA
Plot values of NDVI larger than 0.2
plot(landsatGreen,
main = 'Vegetation cover',
col = "darkgreen",
axes = FALSE, box = FALSE, legend = FALSE)
Create a matrix to be used as an argument in the
reclassify() function:
reclass_m <- matrix(c(-Inf, 0, 1, #water
0, 0.2, 2, #urban
0.2, Inf, 3), #veg
ncol = 3, byrow = TRUE)
reclass_m
## [,1] [,2] [,3]
## [1,] -Inf 0.0 1
## [2,] 0.0 0.2 2
## [3,] 0.2 Inf 3
Classify land cover using the defined threshold values:
landsatCover <- reclassify(landsatNDVI, reclass_m)
Plot the land cover classes:
Save the reclassified raster landsatcover in the GeoTiff
format:
writeRaster(landsatCover,
filename = "clean_data/landsat_landcover.tif",
overwrite = TRUE)
library(landscapemetrics)
library(landscapetools)
#landsatCover <- raster('clean/landsat_landcover.tif') #reload saved raster
show_landscape(landsatCover, discrete = TRUE)
Check if the raster data is in the right format:
check_landscape(landsatCover)
## layer crs units class n_classes OK
## 1 1 projected m integer 3 ✔
E.g. Area of each patch in the landscape:
lsm_p_area(landsatCover)
## # A tibble: 12,551 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 patch 1 1 area 14.6
## 2 1 patch 1 2 area 0.09
## 3 1 patch 1 3 area 27.3
## 4 1 patch 1 4 area 0.09
## 5 1 patch 1 5 area 0.18
## 6 1 patch 1 6 area 0.36
## 7 1 patch 1 7 area 0.09
## 8 1 patch 1 8 area 0.09
## 9 1 patch 1 9 area 0.09
## 10 1 patch 1 10 area 0.09
## # ℹ 12,541 more rows
E.g. For each class, the total area of all patches:
lsm_c_ca(landsatCover)
## # A tibble: 3 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 1 NA ca 6737.
## 2 1 class 2 NA ca 28504.
## 3 1 class 3 NA ca 42930.
E.g. For each class, the average area of patches:
lsm_c_area_mn(landsatCover)
## # A tibble: 3 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 1 NA area_mn 3.23
## 2 1 class 2 NA area_mn 4.52
## 3 1 class 3 NA area_mn 10.3
E.g. Total area of the landscape (all three land cover classes):
lsm_l_ta(landsatCover)
## # A tibble: 1 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA ta 78170.
Spatial data used in this document:
© XP Song